A phenomenological glass model for vibratory granular compaction 
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A model for weakly excited granular media is derived by combining the free volume argument of 
Nowak et al. [Phys. Rev. E 57, 1971 (1998)] and the phenomenological model for supercooled liquids 
of Adam and Gibbs [J. Chem. Phys. 43, 139 (1965)]. This is made possible by relating the granular 
excitation parameter V, defined as the peak acceleration of the driving pulse scaled by gravity, to 
a temperature- like parameter 77(F). The resulting master equation is formally identical to that of 
Bouchaud's trap model for glasses [J. Phys. I 2, 1705 (1992)]. Analytic and simulation results 
are shown to compare favourably with a range of known experimental behaviour. This includes 
the logarithmic densification and power spectrum of fluctuations under constant r\, the annealing 
curve when n is varied cyclically in time, and memory effects observed for a discontinuous shift 
in 77. Finally, we discuss the physical interpretation of the model parameters and suggest further 
experiments for this class of systems. 

PACS numbers: 05.40.-a, 45.70.Cc, 64.70.Pf 
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I. INTRODUCTION 

It is well known that, with the appropriate driving and 
boundary conditions, granular matter can approximate 
each of the three major states of matter: gas, liquid and 
solid ELpJ. Conspicuous by its absence is a glass state; 
that is, a state where the relaxation times far exceed 
the observational timeframe |^-0]. However, it is be- 
coming increasingly clear that the granular analogue of 
glass has been found in a recent series of experiments per- 
formed at the University of Chicago pHll]]- They mea- 
sured the density of a system that was weakly perturbed 
or 'tapped' by the application of a driving pulse to the 
container. A first indication of glass-like relaxation pro- 
cesses came from analysis of the density p(t) , where t is 
the number of times the system had been tapped, which 
was found to increase only logarithmically slowly M , 



p(t) = Pi 



Ap 



l+Sln(l + t/r) 



(1) 



The fitting parameters pf , Ap, B and r are functions of 
the control parameter T, defined as the peak accelera- 
tion of the driving pulse scaled by gravity, T = a max / g. 
Subsequent experiments in which T was varied during a 
run also behaved in a similar manner to glasses under 
a variable temperature f^Hjj]§-[ll| , suggesting a relation- 
ship between T and some elusive temperature-like quan- 
tity. 

Theoretical attempts to understand the experiments 
have ranged from the construction of toy microscopic 
models to higher level, coarse grained descriptions 
fL2| fzsf . The general consensus has been that the slow 
relaxation is due to frustrated dynamics resulting from 
excluded volume effects. More insightful are the free vol- 
ume arguments postulate d by the Chicago group Q and 
Boutreux and de Gennes pfl], which derive the logarith- 



mic compaction with only a small number of assump- 
tions. Provocatively, these assumptions are also key com- 
ponents in established phenomenological models of glass- 
forming liquids, namely those of Adam and Gibbs p7|] 
and Cohen and Turnbull | p8[ , respectively. This fur- 
ther suggests that the analogy with glasses is a valid one. 
However, both of the granular free volume descriptions 
currently lack any mention of the experimental control 
parameter T and hence must be regarded as incomplete. 

In this paper we demonstrate how one of the free vol- 
ume arguments, namely that of the Chicago group, can 
be expanded into a full model that incorporates T. This 
is made possible by postulating a loose analogy between 
r in the granular system and temperature in supercooled 
liquids, and then using this analogy to incorporate ele- 
ments of the Adam and Gibbs theory. The result of this 
process is a master equation for weakly excited granu- 
lar media that is capable of reproducing a wide range 
of known experimental behaviour. The motivations be- 
hind this work are twofold. Firstly, by focusing on only 
a small number of physical mechanisms, the success of 
the model in emulating the experiments indicates that 
the dominant mechanisms may have been correctly iden- 
tified. It is further hoped that this work may help to 
strengthen the relationship between granular matter and 
glasses. This second goal is easily achieved once we show 
that the derived master equation is identical to that of a 
simple glass model due to Bouchaud f29|-|32| . 

This paper is arranged as follows. In Sec. [n] the 
Chicago group's free volume argument is summarised and 
then expanded to a full model by importing elements of 
the Adam and Gibbs theory. The resulting master equa- 
tion that describes the evolution of the system in time is 
specified. Numerical integration of this master equation, 
plus analytical results wherever possible, are compared 
to the experimental data in Sec. Ill, Of particular im- 



portance here is an explanation for the apparent contra- 
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diction between the experiments and any model based 
on the free volume approach, concerning the supposed 
T-dependence of the projected final density, pi in (Q). 
Further discussion on the physical interpretation of the 
model parameters is given in Sec. [V, as well as suggested 
ways in which the various assumptions behind the model 
may be more rigorously checked. Finally, wc summarise 
our findings in Sec. [y] and make some tentative predic- 
tions for future experiments that may help to further 
elucidate the relevant physical mechanisms in granular 
compaction. 



II. DESCRIPTION OF THE MODEL 

The relationship between the Chicago group's free vol- 
ume description and the Adam and Gibbs theory is that 
they both regard the dominant relaxation process to be 
the cooperative rearrangement of particles. The corre- 
spondence between the two theories can be taken further 
by postulating the existence of a temperature-like noise 
parameter rj(T) for weakly excited granular matter. This 
procedure forms the basis of our work, and is described 
in full below. For current purposes it is sufficient to pro- 
vide a somewhat heuristic description of the model; a 
fuller discussion of the various parameters can be found 
in Sec. fvl. 



A. First principles derivation 

The Chicago group's ||, and also Boutreux and de 
Gennes' arguments employ the concept of the mean 
free volume per particle, here denoted Vf . For a system 
of N particles occupying a total volume V, V[ is defined 

as 



_ V-Vmta I 1 

Vf = = v<r \ 

N ° \p Pmax 



(2) 



where v g is the volume of a single particle and p is the 
volume fraction Nv g /V. Units are chosen so that the den- 
sity of a single grain is unity, hence p is also the density 
of the system. Following |2q1 , p max = Nv g /V m in is identi- 
fied with the most compact state possible in a disordered 
system, i.e. the random close-packing limit. In what 
follows we shall fix p max = 0.64, believed to be the ran- 
dom close packed density for a system of monodisperse 
spheres p3[ . 

The Chicago group postulated that the compaction 
process is dominated by the cooperative rearrangement 
of local domains of particles. If z is the number of par- 
ticles in a region that can rearrange independently of its 
environment, they argued there is a lower cut-off 



below which there is not enough free volume available to 
allow reconfiguration. Roughly speaking, z* is the num- 
ber of particles that, by adding up their individual free 
volumes, can make a single 'hole' big enough to allow ex- 
actly one particle to fit through. The explicit dependence 
of z* on p is found by combining @ and (||), 



P Pt> 



(4) 



By assuming that the density increases at a rate propor- 
tional to e~ z , it is now possible to derive the logarithmic 
compaction law p(t) ~ — l/ln(t) fl. 

As mentioned in the introduction, the theory just out- 
lined is incomplete as it does not incorporate the exper- 
imental control parameter T. In an attempt to resolve 
this deficiency, we observe that a similar description for 
cooperative relaxation is also central to the theory pro- 
posed by Adam and Gibbs for structural relaxation in 
supercooled liquids ]27| . An intermediate stage of their 
calculations is of interest here, namely that the relaxation 
rate W can be expressed as a function of temperature T 
as 



W(T) oc exp 



z*AE 



(5) 



where AE is the free energy barrier per particle, fee is 
Boltzmann's constant and z* is again the smallest num- 
ber of particles that can rearrange independently of their 
environment (which was ultimately related to the config- 
urational entropy). 

The principle assumption behind our current work is 
that an expression analogous to (||) also holds for weakly 
excited granular media. More precisely, we propose that 
a region with local density p reconfigures at a rate 



W(p, T) oc exp 



V(T) J 



(6) 



where z* is related to p via (f|). AE can be interpreted as 
a gravitational potential energy barrier per particle, and 
ry(r) gives some measure of the degree of excitation of the 
system. Note that although 77(F) plays the role of fee T, 
we stop short of referring to it as a 'granular tempera- 
ture' and instead regard it as a noise parameter which is 
defined by (||), with the only restriction that 77(F) should 
be a monotonic increasing function of T. In what follows 
77(r) is essentially treated as a fitting parameter. The 
physical meaning of 77(F) and AE is discussed further in 
Sec. 



IV 



z> z* 

V{ 



(3) 



To fully specify the model, some rule is required that 
gives the density of a region after it has reconfigured. In 
general this will depend on its density before reconfigu- 
ration as well as 77(F), but for simplicity we shall ignore 
such considerations here and simply assume that the den- 
sity after reconfiguration is given by the fixed probability 
density function fi'(p). Specifically, p'(p)dp is the prob- 
ability that a region 'falls' into a configuration with a 
density in the range [p, p + dp) . The prior distribution 
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fi'(p) (re-expressed in terms of the total energy barrier 
E - see below) will play a central role in our model, al- 
though it shall be demonstrated that, over timescales rel- 
evant to the experiments, the model is essentially robust 
to the particular choice of p'(p). This is fortuitous, as 
the precise form of fjf(p) is unknown and we have instead 
considered a range of plausible functional forms. 



B. Summary of the model 

Since the reconfiguration rate W(p,T) given in (||) 
depends on p only via the total energy barrier E = 
z*(p)AE, it is convenient to now make the change of 
variables p — > E, where 



E = z*AE = A[ — 

V P /°max 

A = aAE , 



(7a) 
(7b) 



which is one-to-one and hence invertible for all p G 
[0, Pmax) and E £ [0, oo). Thus, in what follows, the state 
of the system at any given time t will in the first place 
be defined by the distribution of energy barriers P(E, t), 
and only then shall the mean density p(t) be found by 
inverting the mapping (M) and averaging over P(E,t), 



P(t) = / 
Jo 



P(E,t) 

A 
E 



+ 



dE 



(8) 



Pn 



Note that, in principle, small values of E should be dis- 
allowed to reflect the fact that low density configurations 
are not mechanically stable and will not arise. For the 
sake of simplicity we choose to ignore this subtlety here. 

The master equation for P(E, t) can be derived as fol- 
lows. The rate at which a region with a local barrier E{p) 
reconfigures is given by ujo eT E ^ 71 ', where the constant uq 
fixes the timescale. After reconfiguring, the region falls 
into a state with a new barrier E new with a probabil- 
ity p(E nevr ), where p(E) is just p!(p) after the change of 
variables, p(E)dE = p'(p)dp. Assuming that the num- 
ber of taps t can be well approximated as a continuous 
variable, P(E, t) evolves in time according to 

i dP(E,t) _ E/v 



LUq dt 



= - e-WPiE, t) + w(t) p(E) , (9a) 



p OO 

,(t) = / c- E ^P(E,t)dE . 
Jo 



(9b) 



The first and second terms on the right hand side of 
( |9a| ) correspond to regions with barriers E before and 
after a reconfiguration event, respectively. Conservation 
of probability is ensured by u>(t), which is the total rate 
of reconfiguration events at time t. 

Remarkably, the coupled equations ( |9a| ) and ( |9~b| ) are 
identical to the trap model of Bouchaud, which is known 
to qualitatively reproduce many features of spin glasses 



and supercooled liquids (29 32 . Thus the model we have 
derived can also be viewed as Bouchaud's trap model, 
with a mapping from the energy barrier E to density p 
that is reached via the two-stage process of first assuming 
that E is proportional to the smallest region that can re- 
arrange independently of its environment, a la Adam and 
Gibbs, and then using the Chicago group's free volume 
argument to relate the size of this region to its density. 
The relationship with Bouchaud's trap model is useful 
as it allows known analytical results to be transferred to 
this application, as described in the following section. 



III. COMPARISON TO THE EXPERIMENTS 

In this section we compare the predictions of the model 
to the experimental results given in [js|-|i"l|]. The general 
procedure employed throughout was to numerically in- 
tegrate P(E, t) in time according to the master equa- 
tion (||) from an initial state P(E,0), using the method 
described in Appendix |X| Ideally P(E, 0) would be cho- 
sen to mimic the distribution of density in the apparatus 
after the preparation phase, but since such information 
is not available we have instead employed the natural 
choice of P(E,0) = p(E), which formally corresponds to 
an instantaneous 'quench' from rj — oo. No significant 
deviations are expected for other initial conditions after 
an initial transient. Once P(E, 0) was fixed, the constant 
A in (0) was chosen by trial-and-error to give an initial 
density close to the experimental value p(0) w 0.58. The 
density p(t) was extracted at regular intervals by numer- 
ical evaluation of (||). 

Each simulation was repeated for two different choices 
of p{E), namely an exponential p(E) — -^e~ E / E ° and a 

Gaussian p(E) = \j2j ■Ko 2 <ir E / 2a , where without loss 
of generality we now choose units such that Eg = a = 1 . 
Other p(E) were also considered for the compaction un- 
der constant r\ described in Sec. Ill A and were found to 
give the same behaviour for t ~ 10 4 taps, indicating that 
the model is robust to the particular choice of fi,(E) over 
the experimental timeframe. However, this robustness 
does not extend to the t — > oo limit, where it is already 
known that different p(E) can give qualitatively different 
behaviour. This is discussed thoroughly in |p2[] , but in 
brief, an exponential tail p(E) ~ e~ E gives rise to a glass 
transition at r] = 1 , in the sense that an equilibrium so- 
lution only exists for 77 > 1 . This can be seen by simply 
setting dP/dt — in the master equation (0), 



Peqm(-B) 



lim P(E, t) 

t— >oo 



Lo{oo)e E/ri p{E) , (10) 



which is not normalisable for rj < 1 if p(E) ~ e~ E , 
and hence equilibrium cannot be reached. By contrast, 
if pi{E) decays more rapidly than exponentially, e.g. 
if it has a Gaussian tail, then an equilibrium solution 
exists for all 77 > 0, although the equilibration time 
may be excessively large for small 77. Note that this 
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model quite generally predicts that the limiting density 
Poo = hm^oo p(t) is a monotonic decreasing function 
of rj. A proof of this is given in Appendix p[ 



A. Constant excitation intensity 

Simulation results for the mean density pit") over a 
range of r\ is given in Fig. [l]. Also given are fits to the 
empirical law (|l]), demonstrating that it is well obeyed 
with either an exponential and Gaussian p(E). We have 
also checked and found similar logarithmic behaviour for 
a selection of other p(E), such as uniform on [E ,Ei], 
both with Eq = and Eq > 0, Gaussian with a non- 
zero mean, Cauchy, and exponential limited to the range 
[Eq, E{\. However, logarithmic relaxation is not expected 
for pathological p(E) such as S(E — Eq) or exp(— e E ). 



°- 0.61 




1000 
number of taps) 



10000 



FIG. 1. Plot of p(t) versus t on log-linear axes for a Gaus- 
sian p(E) = y/2/ne~ E /2 with p max = 0.64, lo = 0.1 and 
A — 0.053. From bottom to top, the lines correspond to 
n = 0.005, 0.03, 0.1, 0.2 and 0.3 respectively. The circles 
are fits to the empirical law ([j]). (Inset) The corresponding 
results for an exponential p(E) = e~ E and A — 0.05, with 
n = 0.002, 0.02, 0.1, 0.2 and 0.5 . 

The logarithmic behaviour can be understood by con- 
sidering the scaling solution to the master equation (||) 
already found by Monthus and Bouchaud [B2[. They 
demonstrated that, after a short transient, P[E, t) can 
be expressed in terms of a single scaling variable u as 



P[E,t) = -u<p(u) 



OJot 



(11) 



Strictly speaking this is only true for an exponential p(E) 
below the glass point, but it was also demonstrated that 
a Gaussian p(E) admits a similar scaling solution until 
a time t* ~ u>Q 1 exp(l/r] 2 ), which may lie well beyond 



the experimental timeframe when rj is small. The phys- 
ical picture underlying this scaling behaviour is that the 
sizes of the cooperatively rearranging regions, which are 
proportional to E = r\ ln(a;niii), are increasing logarith- 
mically in time. A logarithmic increase in domain size 
has also been found in the Tetris model 

Over timescales for which the scaling solution (|Tl| ) 
holds, the density can be expressed in terms of <f>(u) by 
changing variables from E to u in (||), 



P(t) 

Pmax 



= l - 



4>(u) 



=fc i 



Ap n 



■ ln(u)otu) 



■ du 



(12) 



The similarity of this expression to the empirical law (Q) 
is striking. The primary difference is that here we must 
integrate over a distribution of u, which will in general 
introduce corrections to the simple logarithmic law. The 
simulation results in Fig. |l| demonstrate that any such 
corrections are at most small. 

The form of the theoretical prediction (jl^) makes it 
difficult to calculate the fitting parameters Ap, B and t 
in the empirical law (|l|). However, one parameter that 
can trivially be fixed is the projected final density pf , 
which is always equal to p max here, regardless of r/(F). 
In contrast, the experiments seem to indicate that p{ 
is a non-monotonic function of T ||. There is no easy 
way to resolve this discrepancy. For instance, one can- 
not simply assume that /? ma x is itself a function of T, 
i-e. pmax = Pmax(r). Quite apart from the conceptual 
difficulties this would invoke for the physical meaning of 
Pmax ! it would allow situations in which negative free 
volume could arise, for instance by first allowing a sys- 
tem to relax arbitrarily close to p m ax(r) and then sud- 
denly changing to a T 1 for which p ma x(r') < p max (r). 
By definition, Vf would then be negative. Note that this 
contradiction is not specific to this model but will arise 
whenever the definition of free volume (||) is used. 

We believe the solution to this problem lies in the 
range of t over which the data fitting has been performed. 
As mentioned previously, the scaling solution (|ll]), and 
hence the logarithmic relaxation, only applies after a 
short transient, typically t ~ 10 2 — 10 3 taps. However, 
we have found that it is still possible to attain a very 
reasonable fit to the empirical law over the whole range 
< t < 10 4 , but only at the expense of predicting the 
wrong pf . This is clearly demonstrated in Fig. ||[ which 
shows that a fit that works well for < t < 10 4 fails when 
extrapolated to larger t, whereas fixing pf = p max gives 
an initially poorer fit but recovers the correct asymptotic 
behaviour. Transferring this insight to the experiments 
suggests that discarding the first 1-10% of the experimen- 
tal data points and then repeating the fitting procedure 
would result in a similar logarithmic compaction law as 
before, but with pf independent of T. The various time 
regimes in this model are summarised schematically in 
Fig. |. 
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FIG. 2. Comparison between different choices of fitting 
parameters over different ranges of t. The thick line 
is the density p(t) for a Gaussian fi(E) at 77 = 0.1 
with p max = 0.64 . The circles correspond to the fit 
p(t) = 0.63-0.053/(1 + 2.3 ln(l+t/16)) and the crosses corre- 
spond to p(t) = p max — 0.06/(1 + 0.241nt). (Inset) The same 
plots over the experimental timeframe 1 < t < 10 4 . 
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FIG. 3. Schematic of densification under constant 77. The 
line becomes straight, indicating that the density is loga- 
rithmically increasing towards p max , once the system en- 
ters into the scaling regime, which overlaps with the experi- 
mental window. The scaling regime continues indefinitely if 
p(E) ~ e~ E ^ Ve and 7] < 1]g. For 77 > % , or for a p(E) with 
a tail that decays faster than exponentially, the scaling be- 
haviour ceases at some late time and the density reaches its 
equilibrium density p x < p max • 



sity to vary in time, which roughly corresponds to varying 
the temperature in other slowly relaxing systems 
Two time-dependencies will be considered in this paper. 
The first is the 'annealing curve,' which was experimen- 
tally attained by cyclically ramping T in a stepwise fash- 
ion between some high value T = T\ and T = [p|,^0[. 
Slowly decreasing T removes low density local configura- 
tions without creating many new ones, hence the term 
'annealing.' The second protocol for varying V will be 
investigated in the next subsection. 

The annealing curve for this model is obtained by al- 
lowing 77 to smoothly vary from to some value 771 to 
to 771 again, where the duration of each leg is denoted 
by tics ■ Simulation results for 7ji cg = 10 6 are given in 
Fig. ||. The experimental annealing curve has a sim- 
ilar shape, except that the initial density increase for 
small r is noticeably slower than that for small 77 ||. 
This may simply be due to a non-trivial mapping from 
T to 77, as discussed in Sec. [V. Note that the second 
and third legs in Fig. ^ form a reversible curve which is 
nonetheless out of equilibrium for small 77. Observe also 
the presence of a narrow hysteresis loop, which is also 
present in microscopic models Jl4],[l7],^2| but has never 
been systematically searched for in experiments. The 
area of this hysteresis loop decreases for slower cooling 
rates, as demonstrated in Fig. @. 
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FIG. 4. The annealing curve for a Gaussian p(E). The 
symbols refer to the initial increase in 77 up to 77 = 0.6 (open 
squares), the decrease to 77 = (filled circles) and the second 
increase (crosses). For each leg, 77 was varied smoothly over 
hag = 10 6 taps. The thick line is the equilibrium density. (In- 
set) The same for an exponential p(E) over a range of 77 that 
includes the glass point. 




B. Annealing curve 

Further insight into the nature of the system's relax- 
ation properties can be gained by allowing the tap inten- 



To interpret these results in a glassy context, re- 
call that the initial conditions were chosen to conform 
to the equilibrium state at 77 = 00, i.e. P(E,0) = 
-Peqm(-E)|r;=oo = l^(E). This would be valid if the ini- 
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tial low density configuration in the experiments corre- 
sponded to an equilibrium state for very large tapping 
intensity T, which seems plausible. Thus the start of the 
first leg corresponds to a rapid quench from high 77 to 
X] « 0, leaving the system far from equilibrium. The rate 
of compaction is initially rapid but slows as the density, 
and hence the relaxation times, increase. For sufficiently 
high 77, the density reaches and starts to follow the equi- 
librium curve, rapidly erasing the memory of its history. 
As 77 is lowered a second time, this time corresponding 
to a slow quench, the system remains near equilibrium 
until some value r\§ (which depends on the cooling rate) 
when the relaxation time rapidly increases and the sys- 
tem essentially freezes. Thus the difference between the 
first and third legs in Fig. || can be understood as the 
recovery from a rapid and slow quench, respectively. 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 
r| (arb. units) 



FIG. 5. Variation of the annealing curves with cooling rate 
for a Gaussian fi(E). Results are given for a total time per 
leg of ti e g = 10 7 (thin solid line), 10 5 taps (dashed line) and 
10 3 taps (dot-dashed line). The second and third legs are 
reversible for ti eg ~ 10 5 . The thick line is the equilibrium 
density. 



C. Shift in the excitation intensity 

Recent experiments have investigated the effect of al- 
lowing r to 'shift' from a constant value Tq to another 
constant value Ti at a given time to []lT| . It was found 
that the system evolved in a way that depended on its 
history as well as its current density and excitation in- 
tensity r, representing a form of memory akin to that 
in glassy systems J|||. Plotted in Fig. g are the corre- 
sponding results for this model, where 77 = 770 = 0.5 until 
t = to = 50 taps, when it changes to rjx = 770 + A77. The 
sign of the initial density change is opposite to the sign of 
A77, as in the experiments, although this is not entirely 
general and the behaviour is reversed if to is too small. 
Also shown in the inset is the case when 77 is changed 



from different 770 to the same value 771 = 0.3 when the 
density reaches a predetermined value. Again there is 
qualitative agreement with the experiments. 

The analogy with glass systems suggests that the 
timescale of the response to a shift in 77 at to should 
scale with to in some manner |34j]. With this insight, we 
now make the following prediction, in the hope it may be 
tested experimentally. Let Ap(t — to) be the difference in 
density at time t between the perturbed system and an 
unperturbed one, i.e. one with A77 = . Plotted in Fig. ^ 
is Ap(t — to) for a shift from a low to a high 77 at times 
to = 10 3 , 10 4 , 10 5 , 10 6 and 10 7 . In each case there is a 
well-defined time for the peak response t resp , which in- 
creases with to . Known results for the trap model with 
an exponential p{E) suggest that t rcsp ~ t Q l0 ^ Vl , where 
the exponent is independent of to |55| . We find this to be 
a good first approximation to our data, as demonstrated 
by the inset to Fig. |?], although there are corrections aris- 
ing from the non-linear mapping from E to p, which dis- 
torts the underlying scaling behaviour in E. There are 
also additional small corrections when using a Gaussian 
fx(E). 

Finally, the experiments briefly investigated what hap- 
pens when r is allowed to return to its initial value after 
St taps at a higher value Ti pH^ . For comparison, the 
equivalent results from this model are given in Fig. ||. 
The observed trend is in accord with the experimental 
observations. A full study of this variation in ?7(t) for all 
Vo 1 to and St is beyond the scope of this paper and 
will not be discussed further here. 
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FIG. 6. Response to a shift from 77 = 0.5 to 77 = 0.5 + A77 
at a time t — 50 for a Gaussian p.(E). From top to bottom, 
the lines refer to A77 = —0.3, -0.2 and -0.1 (dashed lines), 
(thick line) and 0.1, 0.2 and 0.5 (thin solid lines). (Inset) Here 
77 = 0.1 (solid line), 0.3 (dashed line) or 0.5 (dot-dashed line) 
until the first time to when p(to) > 0.61, afterwhich 77 is fixed 
at 0.3 in each case. 
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FIG. 7. Response to a shift from r\ = 0.1 to r\ = 0.5 at 
a time to for a Gaussian p(E), where Ap(t — to) is the den- 
sity difference between perturbed and unperturbed systems. 
From top the bottom on the right hand side, the lines refer 
to t = 10 3 , 10 4 , 10 5 , 10 6 and 10 7 respectively. An exponen- 
tial p(E) behaves similarly. (Inset) The same data plotted 
against (t — to)Ao , where a = 0.1/0.5 is the ratio of rj before 
and after the shift. 
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FIG. 8. Plot of the recovery from a short time at a higher 
rj for a Gaussian p(E). Here the system has been relaxed for 
a time to = 10 4 taps at rj — 0.1, then held at rj = 0.3 until 
ti = to + St when it reverts back to rj = 0.1 again. From 
bottom to top, the lines refer to St — 1, 2, 4, 8 and 16, re- 
spectively. (Inset) The raw data for St = 16 (dashed line) 
compared to the unperturbed system (solid line). 



D. Fluctuations and power spectra 

In a finite system the density in equilibrium is not con- 
stant but fluctuates about its mean value. To investigate 
density fluctuations in this model, a different version of 
the code was employed which explicitly simulates a sys- 
tem consisting of N separate subsystems (details given in 
Appendix [A]) . Fig. ^ shows the probability distribution 
Q(Ap) of fluctuations Ap = p(t) — p mo dai for N = 500 
and different rj. To first approximation Q(Ap) is Gaus- 
sian, but it is slightly skewed towards lower densities, 
becoming more so as rj is lowered. The skewness arises 
from the non-linear mapping from E to p, which exag- 
gerates fluctuations to lower densities whilst suppressing 
those to high densities. There is also a cut-off for very 
large | Ap | , when P(E, t) has deviated significantly from 
Peqm{E). The experiments exhibited Gaussian fluctua- 
tions with some anomalous deviations for Ap > 
this is discussed below. 




X = (Ap)" sgn(Ap) (arb. units) 



FIG. 9. Fluctuations around the modal density p m odai 
in equilibrium for an N = 500 element system with a 
Gaussian p(E). Q(Ap), the probability of a fluctuation 
Ap = p(i) — p m odai , is plotted against X — (Ap) 2 sgn(Ap) 
on a log-linear plot, which would give a symmetrical triangle 
if Q(Ap) was Gaussian. The different lines refer to rj = 1 
(solid), rj — 0.8 (dashed) and rj = 0.6 (dot-dashed). (Inset) 
The same for an exponential p(E) with rj = 2 (solid line) and 
rj — 1.6 (dashed line). 

The power spectra S(f) of density fluctuations in equi- 
librium for various rj and a Gaussian p(E) is given in 
Fig. [l(| S(f) ~ l// 2 for / greater than a high frequency 
shoulder /h , where /h is only weakly dependent on ry, 
although it should be stressed that this model focuses on 
cooperative relaxation modes and is not intended to de- 
scribe the high frequency, single particle dynamics. For 
low frequencies, S(f) appears to obey non-trivial power 
law behaviour S(f) ~ l//' 5 , where the exponent 6 can 
be very approximately fitted to 6 w 1 — 77 over the range 
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1CP 5 < f < 10~ 3 . However, the analysis given in Ap- 
pendix |C| shows that this is not the true asymptotic be- 
haviour and S(f) — > 1//° as / — > 0. The crossover to 
1//° behaviour occurs around a low frequency shoulder 
fi, , where /l — * rapidly as 77 — ► 0. For an exponential 
there is only one shoulder frequency separating the 
high frequency, l// 2 regime from a low frequency regime 
in which S(f) ~ l//' 5 , where £ = 2 — 77forl<77<2 and 
5 = for 77 > 2. Note that there is no 1//° region for 
77 < 2, even though the system is in equilibrium. This 
apparent anomaly is explained in Appendix |C| 



In the experiments, the power spectra were found to 
obey non-trivial power law behaviour S(f) ~ l//" 5 , with 
S = 0.9 ± 0.2, between two corner frequencies /l and /h 
that both decreased as V was lowered ||. More complex 
behaviour was observed for larger T towards the bottom 
of the apparatus. The results from this model are in par- 
tial agreement; for instance, it is still one of few that can 
exhibit a l/f s regime with S s» 1 (sec also |2~l]| ). There 
are some discrepancies, but these may simply be due to 
processes not currently incorporated into the model, such 
as single particle dynamics or the existence of metastable, 
high-density crystalline domains. 
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FIG. 10. S(f), the power spectrum of frequency /, for 
an JV = 500 system and a Gaussian n(E), at 77 = 0.6 (solid 
line), 77 = 0.4 (dashed line) and 77 = 0.2 (dot-dashed line, also 
vertically shifted by a factor of 500 for clarity). Each S(f) 
was calculated over « 10 8 points. The slopes of the thick 
line segments are indicated. To speed convergence, the ini- 
tial configuration was chosen as the known equilibrium state 
Pc qm (S) for N = 00 ©, although the 77 = 0.2 system was 
still evolving towards its slightly different finite N equilibrium 
state during the run. This accounts for the anomalous steep 
slope for small /. 
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FIG. 11. S(f) in equilibrium for an exponential fi(E) and 
N = 500, with 77 = 1.7 (solid line, also shifted vertically by 
a factor of 10) and 77 = 1.3 (dashed line). The analytical 
predictions from Appendix |c] are indicated by the thick line 
segments. 



IV. DISCUSSION OF THE MODEL 
PARAMETERS 



Given the success of this model in reproducing the 
experimental phenomenology, it is natural to ask if its 
principle assumptions can be placed on a firmer founda- 
tion. Inparticular, a number of parameters introduced in 
Sec. II A have so far been treated somewhat heuristically. 
To redress the balance, we now discuss the physical inter- 
pretation of some of these parameters. A more thorough 
analysis may be possible by detailed comparison with a 
microscopic model, for instance. 



1. The noise parameter 77(F) 



It was stressed during the derivation of this model that 
the noise parameter 77 need not bear any relation to the 
concept of granular temperature |36||3^]. By the same 
token, the use of the term 'equilibrium' to describe the 
statistical steady state merely refers to a dynamic equi- 
librium, without supposing any analogy with a thermody- 
namic one. Instead, 77 was defined in the broadest sense 
of simply giving some measure of the degree of excitation 
of the system during a single tap. This loose definition 
makes finding the precise relationship with T difficult. 
Nonetheless it is still possible to predict the overall shape 
of 7y(r), as we now argue. 

For 77 to be non-zero, the particles must at the very 
least separate from their nearest neighbours. Experi- 
ments on vibrated granular systems often claim to find 
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some critical T c such that the relative motion of the par- 
ticles is either minimal or non-existent for T < T c (usu- 
ally 1 ~ T c ~ 2, see e.g. [p]Pj3^] ). A facile explanation 
for this is to suppose that a granular body is held to- 
gether by frictional forces, and that relative motion be- 
tween adjacent particles is not possible until some static 
friction threshold has been overcome. Since all the nor- 
mal contact forces are proportional to g, the distribution 
of threshold forces will also scale with g and thus the rel- 
evant parameter would indeed be T = a max / g. However, 
friction is not the only relevant mechanism. For instance, 
particle separation will still occur in vertical one dimen- 
sional columns, where friction clearly plays no part. Even 
in this case theory suggests that the relevant parameter 
is again T, at least in the limit of hard spheres |39f| . 

Assuming that a well-defined T c exists, the overall 
shape of rj(T) will obey t](T) sa for small T, only signif- 
icantly deviating from zero for T ~ T c . Just this qualita- 
tive behaviour has been found in simulations of horizon- 
tally vibrated systems |4Cf|. As a corollary, the annealing 
curves presented in Figs. |4| and ^| will be flatter for small 
r when plotted against T rather than rj, in better agree- 
ment with the experimental graphs |9pC[|. 

2. The prior distribution fi(E) 

The Gaussian and exponential p(E) employed in the 
simulations were chosen as plausible first guesses of the 
real [i(E). To calculate the actual n(E) is a non-trivial 
problem, but a first step might be to re-express pi(E) in 
terms of ip(vf), the distribution of free volume after re- 
configuration. Using E = z*AE — Av g /iJf (||), this gives 

, , Ay,, ( Av„ \ 

m = ^i>{-f) ■ as) 

In principle, ip(vf) could be found from a microscopic 
model, such as the parking lot model |P,|l7|-|l9| , for which 
the free volume is also the void volume. 

Note that, from (|l3|), the tail of /i(-E), which is so 
important to the long-time relaxational properties of 
Bouchaud's trap model, can be related to the Vf — > + 
behaviour of ij)(vf). For example, if ip{vi) vanishes ac- 
cording to ip(vf) ~ exp(— a/iJf ), then p(E) will have an 
exponential tail and the trap model predicts a glass tran- 
sition at a finite noise intensity r\ = Av g /a. Similarly, 
ip(v{) ~ exp(— a/vf 2 ) corresponds to a fi(E) with a Gaus- 
sian tail. 



3. The constant A = aAE 

Even though A has been treated as an arbitrary con- 
stant and fixed by the initial conditions, its component 
factors a and AE have a physical interpretation, as we 
now discuss. In (||), a is defined as the constant of pro- 
portionality between z*, the smallest number of particles 



that can cooperatively rearrange, and the ratio v g /vf . In 
the Chicago group's original argument ||, a was set to 1; 
however, we prefer not to fix a at any particular value 
and suggest that it may depend upon particle properties 
such as their shape. For instance, highly irregular par- 
ticles will obstruct motion more effectively than rounder 
particles of the same volume, and so should have a higher 
value of a. 

AE was originally defined as a gravitational poten- 
tial energy barrier. Indeed, assuming that the particles 
interact via hard core repulsion, this is the only avail- 
able potential energy scale in the system. This suggests 
that AE is proportional to the mean vertical displace- 
ment between adjacent particles. If so, then our implicit 
assumption that AE is independent of 77 and p is com- 
patible with the hard sphere Monte Carlo simulations of 
Barker and Mehta jnj, who demonstrated that the dis- 
tribution of contact angles between particles is roughly 
constant over a wide range of shaking amplitudes. 

More importantly, if AE is gravitational potential en- 
ergy, then inspection of the reconfiguration rate (|^) indi- 
cates that rj must also have units of energy, with an en- 
ergy scale that is presumably coupled to the driving. For 
definiteness, suppose 77 can be written as rj — mgAof(T), 
where m is the typical mass of the particles, Aq is the 
amplitude of the driving and f(T) is some function of 
the dimensionless parameter T — a max /g (possibly with 
a threshold around T ~ T c as discussed earlier). Simi- 
larly, write AE ~ mgr 7 where r is the typical particle 
radius. Since AE and rj only appear in the ratio AE/rj, 
m and g will cancel and the dynamics of the model will 
depend on two dimensionless quantities, namely T and a 
dimensionless displacement Ao/r. The existence of a sec- 
ond relevant dimensionless parameter implies that the 
behaviour in response to high amplitude, low frequency 
driving may be qualitatively different from a low ampli- 
tude, high frequency driving with the same value of T. 
This possibility has not yet been explored in the experi- 
ments, which seem to have focused on the low frequency 
regime. 

Note that we could equally have expressed rj in terms 
of the kinetic energy supplied by the driving, i.e. i] = 
towq /(r), where vq is the typical driving velocity. How- 
ever, this is not an independent energy scale as vq can be 
dimensionally related to a max and Aq by vq ~ \J a max Aq . 
We only mention this latter alternative because simula- 
tions often show scaling plots in terms of T and vq (see 
e.g. E3] and references therein). 



V. SUMMARY AND CONCLUSIONS 

To summarise, we have constructed a simple model for 
weakly excited granular media that combines the Chicago 
group's free volume argument with elements of the super- 
cooled liquid theory of Adam and Gibbs. Integration of 
the master equation has shown that the model behaves 
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in a similar manner to the experiments for each of the 
situations considered. Some slight discrepancies remain 
with the power spectra, but these may be due to mecha- 
nisms currently lacking from the model, such as ordering 
effects and crystallinity, depth dependency or wall effects. 
It would be interesting to see if any of these mechanisms 
could be incorporated into an extended version of model. 
It may also be possible to introduce orientational degrees 
of freedom and compare the results to recent experiments 
on nylon rods [j42[ . 

The model has also been used to predict the man- 
ner in which the time of the peak response to a shift 
in r at t — to scales with to , as discussed in Sec. Ill C . 
This prediction could be tested experimentally and may 
help to differentiate between the large number of models 
that have so far been proposed ]l2|-p5|], as it seems un- 
likely that they will all give the same scaling behaviour. 
Further insight into the physical mechanisms underlying 
the compaction process could be gained by measuring 
the typical size of reconfiguring regions as a function of 
time, or by seeing if the locations of such regions are 
spatiotemporally correlated. Such measurements could 
be performed in simulations, or by direct visualisation of 
two dimensional experiments f43|| , for instance. 

Finally, we note that the relationship between granular 
media and glasses can be given a more intuitive appeal 
by the following simple argument. Consider sand poured 
from a great height into a container. When the particles 
first hit the surface of the forming sandpile, the direction 
in which they bounce will essentially be random, giving 
rise to a large random velocity component. This corre- 
sponds to a highly excited state with (in our notation) 
a high 77. However, the particles will rapidly lose their 
kinetic energy by inelastic collisions and will soon come 
to rest, jamming under gravity into a static, disordered 
configuration with 77 = 0. It is not difficult to see how this 
sequence of events can be related to the rapid 'quench' 
of a supercooled liquid or other glass-forming material. 

Just after the initial submission of this work, we be- 
came aware of a master equation for the glass transition 
due to Dyre [Q, which is similar to Bouchaud's equa- 
tion studied in this paper but with a built-in cut-off in 
the range of allowed energies. Also, it has been brought 
to our attenti on th at the two regimes of vibration men- 
tioned in Sec. IV 3 have previously been discussed in the 
context of size segregation by Mehta and Barker Ha] . 
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APPENDIX A: SIMULATION DETAILS 

The bulk of the simulation results were obtained by 
numerical integration of the continuous master equa- 
tion (||). P(E, t) was defined on a mesh of points Py = 
P(iSE,j 6t), where < i < i max and j > 0. Care was 
taken to ensure that E max = imax 

SE was set sufficiently 
high that there was no significant cut-off to P(E,t) for 
large E. To iterate over a single time step St, oj(t) was 
found from numerical integration of (9b) and then as- 
sumed to remain constant over the required time inter- 
val. This allowed the time evolution equation rt9a]) to be 



solved and Pij+i found from Py Vi. The whole distri- 
bution was then renormalised by a factor Pij) to 
correct for the non-conservation of probability resulting 
from the assumption of a constant w(t). For relaxation 
under constant 77, simulation times were improved by em- 
ploying a geometric mesh with a linearly increasing time 
step St oc t. This allowed for times up to t = 10 10 to be 
reached with only modest CPU time. 



For the density fluctuations investigated in Sec. [II D, 



the continuous master equation was of no use and an 
alternative method was employed which explicitly in- 
cluded finite size effects. This involved assigning N ar- 
ray elements a barrier E\ , i — 1 . . . N, according to the 
chosen initial conditions. At every time step St = 1, 
each element was assigned a new barrier with probability 
u>o e~ Ei ^ , where the new barrier values were drawn from 
the prior p{E). The density pi of each element was found 
by inverting the mapping (|t]) , and the mean density cal- 
culated by straightforward summation, p(t) = jj^Pi ■ 



APPENDIX B: MONOTONICITY OF Poo (rj) ON 77 

In this appendix it is shown that the asymptotic den- 
sity Poo = lfm^oo p(t) is a monotonic decreasing function 
of 77 for essentially any fi(E). If an equilibrium state ex- 
ists, it takes the form P eqm (E) = LJooiv) e E ' r, p,(E) and 
hence from (||), 

Woo(r?)e £ /V(£) 



Poo 



1 + E/Ap n 



■dE 



^oo(?7) = nm u(rj,t) 

t — *oo 



e £, 'V(£') dE' 



(Bl) 



(B2) 



Differentiating (Bl) with respect to 77 and rearranging 
gives 



if 



dpooiv) 



_, p(E)p(E')e( E + E 'V r > 



-'max 
r 00 poo 



E/A Pl) 



dEdE' . (B3) 
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After the change of variables u = E + E' and v = E — E' 
and substituting v — > —v over the domain v < 0, the 
right hand side of (jB3|) transforms to 



OO pU 



v 2 e u/ > 



u=0 Jv=0 



U + V 



u — V 



(2Ap max + u + v)(2Ap max +u-v) 



du dv 



(B4) 



Since all the factors inside the integral (B4) are posi- 
tive, it can be trivially deduced that 



dp„ 



drj 



< for all rj. 



(B5) 



Equality is attained in only two cases. The first is if the 
integrand is strictly zero over the entire range v > 0, 
which can only happen in the trivial case of a single val- 
ued distribution p,(E) — S(E — E ). The second and 
more important situation is over a range of i] for which 
no equilibrium state exists. In this case, all the moments 



of P{E, t) diverge as t — > oo and pit) 



fix 



(0). 



A plot of Pooiv) versus r\ for an exponential and a Gaus- 
sian p,(E) is given in Fig. Il2| 



0.64 



0.62 



0.60 



0.58 



3 4 
(arb. units) 



near the glass point has already been derived in |30| ; here 
we consider the / — > + limit in equilibrium for general 
p(E) over a wider range of rj. 

Since each local region is assumed to relax indepen- 
dently of its environment, the total power spectrum S(f) 
is just the spectrum for a single region with a relaxation 
time r averaged over $(r), the distribution of relaxation 
times in equilibrium, 



S(f)<x 



•*(t) 



1 + (2k fr) 



■dr 



(CI) 



The small / behaviour of flCip depends on the asymp- 
totic behaviour of $(t) for large r. If $(r) decays faster 
than r~ 2 , then the / = limit exists and S(f) exhibits 
the expected 1//° noise for low frequencies. However, if 
$(r) - t- x with 1 < x < 2, then S(f) - Uf 2 ~ x , as can 



be readily seen by substituting for fr in (CI) (note that 
x > 1 since $(t) is normalisable) . 

For the trap model, $(t) can be found for any 
given p{E) by simply making the change of variables 
t = ^e B / ?) into the expression for P c(pn (E), equa- 
tion (p"o|). Thus, an exponential pb(E) ~ e~ E ^ ns gives 
$(r) ~ T -"/^, implying that S(f) - 1/ '/ 2 -"/% for 
?y g < 77 < 2?7 g . This confirms that 7^ 1//° for this 
range of 77, even though the system is in equilibrium. The 
usual 1//° behaviour is recovered when 77 > 2?7 g , which 
also applies for all r\ when p,[E) decays faster than ex- 
ponentially. In particular, a Gaussian p(E) ~ e - E2 / 2 < j2 
leads to an equilibrium distribution of relaxation times 
of the form 



$(r) 



■ Iii(wot) 



(C2) 



which is suggestive of a power law with a slowly varying 

2 

exponent x — ln(wo T )- Thus one would expect S(f) 
to exhibit approximate power law behaviour over a wide 
range of /, reverting to 1//° only for frequencies compa- 
rable to the 'largest' relaxation time ojqt* ~ e a / r ' . Any 
attempt to fit S(f) to a power law will give an exponent 
that depends on the range of / considered as well as the 
ratio 77/cr. 



FIG. 12. Plot of poo = limt^oo p(t) against rj for an 
exponential p(E) = e~ E (solid line) and the Gaussian 
H(E) = y/2/^e- E2/2 (dotted line) for p max = 0.64. Note 
that the plateau for rj < 1 in the exponential case corresponds 
to the out-of-equilibrium situation where all the moments of 
P(E,i) diverge as t — > 00. 



APPENDIX C: ANALYSIS OF THE POWER 
SPECTRA 

The small / behaviour of the power spectra of density 
fluctuations S(f) is analytically derived in this appendix, 
which extends th e rang e of the numerical observations 
discussed in Sec. HID. The time-dependent spectrum 
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